function [ibeta] = beta_i(a, b, x)
% incomplete beta function Ix(a,b)
% requires betacf and gammaln functions
%


ibeta = 0;
if(x < 0 | x > 1)
   disp('Bad x in routine betai')
   return
end

if(x == 0 | x == 1)
   bt = 0;
else
   bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
end
if(x < (a+1)/(a+b+2.0))
   ibeta = bt*betacf(a,b,x)/a;
else
   ibeta = 1.0-bt*betacf(b,a,1.0-x)/b;
end
